/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package stendcontrol;

import java.util.ArrayList;

/**
 *
 * @author Asko
 */
public class Integrals {


    public static double findOrdinatesNum(int modelNum, double lambda,
            double omega) {
        double ordNum = 0;
        switch (modelNum) {
            case 1: {
                ordNum = 9;
            }
            break;
            case 2: {
                ordNum = 13;
            }
            break;
            case 3: {
                ordNum = 19;
            }
            break;
            case 4: {
                ordNum = (47.0 * omega) / (2.0 * Math.PI * lambda);
            }
            break;
            case 5: {
                ordNum = (47.0 * omega) / (2.0 * Math.PI * lambda);
            }
            break;
            case 6: {
                ordNum = (47.0 * omega) / (2.0 * Math.PI * lambda);
            }
            break;
            default:
                break;
        }
        return Math.round(ordNum);
    }

    public static double findGamma(double k, double sampIntrv) {
        return 0.4 / ((2 * k + 1) * sampIntrv);
    }

    public static double findInverseLegNorm(int k, double gamma) {
        return 1 / 2 * gamma * (2 * k + 1);
    }

    public static double findSamplingInterval(double omega,
            double error, int modelNum, double lambda) {
        double derivative2 = 0;
        switch (modelNum) {
            case 2: {
                derivative2 = -lambda * lambda;
            }
            break;
            case 3: {
                derivative2 = 3 * lambda * lambda;
            }
            break;
            case 4: {
                derivative2 = lambda * lambda - omega * omega;
            }
            break;
            case 5: {
                derivative2 = -lambda * lambda - omega * omega;
            }
            break;
            case 6: {
                derivative2 = -lambda * lambda - omega * omega;
            }
            break;
            default:
                break;
        }
        return Math.sqrt(8 * error / Math.abs(derivative2));
    }

    public static double lager(double t, double gm, double k) {
        double res = 0;
        if (k == 0) {
            res = Math.exp(-gm / 2 * t);
        } else if (k == 1) {
            res = Math.exp(-gm / 2 * t) * (1 - gm * t);
        } else {
            for (int i = 0; i < k + 1; i++) {
                res += (factorial(k) * Math.pow(-gm * t, i)
                        * Math.exp(-(gm * t / 2.0)))
                        / (factorial(i) * factorial(k - i)
                        * factorial(i));
            }
        }
        return res;
    }

    public static double lager0(double t, double gm) {
        return Math.exp(-gm / 2 * t);
    }

    public static double lager1(double t, double gm) {
        return Math.exp(-gm / 2 * t) * (1 - gm * t);
    }

    public static double lagerI0(double t, double gm) {
        return -2 / gm * Math.exp(-gm / 2 * t);
    }

    public static double lagerI1(double t, double gm) {
        return Math.exp(-gm / 2 * t) * (2 * t + 2 / gm);
    }

    public static double lagerIT0(double t, double gm) {
        return (-4 * Math.exp(-gm / 2 * t) * (gm * t / 2 + 1)) / (gm * gm);
    }

    public static double lagerIT1(double t, double gm) {
        return Math.exp(-gm / 2 * t) * (6 * t / gm + 12
                / (gm * gm) + 2 * t * t);

    }

    public static double recLagI(double t, double gm, double k) {
        double sum = 0;
        if (k == 0) {
            return lagerI0(t, gm);
        } else if (k == 1) {
            return lagerI1(t, gm);
        } else {
            for (int i = 0; i < k; i++) {
                sum += recLagI(t, gm, i);
            }
            return -2 * sum - 2 * lager(t, gm, k) / gm;
        }
    }

    public static double recLagIT(double t, double gm, double k) {
        if (k == 0) {
            return lagerIT0(t, gm);
        } else if (k == 1) {
            return lagerIT1(t, gm);
        } else {
            return -((k + 1) / gm) * recLagI(t, gm, k + 1)
                    + (2 * k + 1) * recLagI(t, gm, k) / gm - k
                    * recLagI(t, gm, k - 1) / gm;
        }
    }

    public static double getRecIntegral(double t, double gm, double k,
            boolean leg, boolean withT) {
        if (leg) {
            if (withT) {
                // System.out.println("Leg with t");
                return legRecIntegralT(t, gm, k);

            } else {
                // System.out.println("Leg without t");
                return legRecIntegral(t, gm, k);
            }
        } else {
            if (withT) {
                // System.out.println("Lag with t");
                return recLagIT(t, gm, k);
            } else {
                //  System.out.println("Lag without t");
                return recLagI(t, gm, k);
            }
        }
    }

    public static double leg(double t, double gm, double k) {
        double res = 0;
        if (k == 0) {
            res = Math.exp(-t * gm);
        } else {
            for (int i = 0; i < k + 1; i++) {
                res = res + factorial(k) / (factorial(k - i) * factorial(i))
                        * factorial(k + i) / (factorial(k) * factorial(i))
                        * Math.pow(-1, i) * Math.exp(-(2 * i + 1) * gm * t);
            }
        }
        return res;
    }

    //Интеграл 0 порядка без т  
    public static double legI0(double t, double gm) {
        return -1 / gm * Math.exp(-gm * t);
    }
    //    return  2*Math.exp(-3 * gm * t)/3*gm - Math.exp(- gm * t)/gm ; 
    //Интеграл 1 порядка без т

    public static double legI1(double t, double gm) {
        return -1 / gm * (1 - 2.0 / 3.0 * Math.exp(-2 * gm * t)) * Math.exp(-gm * t);
    }
    //Интеграл 0 порядка с участием т

    public static double legIT0(double t, double gm) {
        return -(t / gm + 1 / (gm * gm)) * Math.exp(-gm * t);
    }
    //Интеграл 1 порядка с участием т 

    public static double legIT1(double t, double gm) {
        return -(t / gm + 1 / (gm * gm)) * Math.exp(-gm * t)
                + 2 * (t / (3.0 * gm) + 1.0 / (9.0 * gm * gm))
                * Math.exp(-3.0 * gm * t);
    }
    //Рекурсивный поиск интеграла без т

    public static double legRecIntegral(double t, double gm, double k) {
        if (k == 0) {
            return legI0(t, gm);
        } else if (k == 1) {
            return legI1(t, gm);
        } else {
            return 2 * (2 * k - 1) / (2 * k + 1)
                    * legRecIntegral(t, gm, k - 1)
                    - (2 * k - 3) / (2 * k + 1)
                    * legRecIntegral(t, gm, k - 2)
                    - 1 / (gm * (2 * k + 1)) * (leg(t, gm, k)
                    - leg(t, gm, k - 2));
        }
    }
    //Рекурсивный поиск интеграла с участием т

    public static double legRecIntegralT(double t, double gm, double k) {
        if (k == 0) {
            return legIT0(t, gm);
        } else if (k == 1) {
            //   System.out.println("In >>>>");
            return legIT1(t, gm);
        } else {
            return 2 * (2 * k - 1) / (2 * k + 1)
                    * legRecIntegralT(t, gm, k - 1)
                    - (2 * k - 3) / (2 * k + 1)
                    * legRecIntegralT(t, gm, k - 2)
                    - 1 / (gm * (2 * k + 1))
                    * (t * (leg(t, gm, k) - leg(t, gm, k - 2))
                    - legRecIntegral(t, gm, k)
                    + legRecIntegral(t, gm, k - 2));
        }
    }

    public static double[] findA(double p[], double b[], double t[]) {
        double a[] = new double[b.length];
        for (int i = 0; i < a.length; i++) {
            a[i] = p[i] - b[i] * t[i];
            //       System.out.println("a[" + i +  "] = " +a[i] );  
        }
        return a;
    }
    //Находим Bi

    public static double[] findB(double[] p, double sampIntrv) {

        double b[] = new double[p.length - 1];

        for (int i = 0; i < p.length - 1; i++) {

            b[i] = (p[i + 1] - p[i]) / sampIntrv;
            //       System.out.println("b[" + i +  "] = " +b[i] );   
        }

        return b;

    }
    //Находим массив т-итых

    public static double[] findT(double numOrd, double sampIntrv) {

        double t[] = new double[(int) numOrd];
        double iter = sampIntrv;

        for (int i = 0; i < t.length; i++) {

            t[i] = sampIntrv - iter;

            sampIntrv += iter;
            //      System.out.println("T[" + i + "] = "  + t[i]);
        }
        return t;

    }
    //Поиск бетта-коффициентов

    public static double[] findBetta(double t[], double a[],
            double b[], double k, double gm, boolean isLeg) {
        double[] Betta = new double[(int) (k)];
        double sum = 0;
        for (int i = 0; i < k; i++) {
            for (int j = 0; j < t.length - 1; j++) {
                if (isLeg) {
                    sum += a[j] * (getRecIntegral(t[j + 1], gm, i, true, false)
                            - getRecIntegral(t[j], gm, i, true, false))
                            + b[j] * (getRecIntegral(t[j + 1], gm, i, true, true)
                            - getRecIntegral(t[j], gm, i, true, true));
                } else {
                    sum += a[j] * (getRecIntegral(t[j + 1], gm, i, false, false)
                            - getRecIntegral(t[j], gm, i, false, false))
                            + b[j] * (getRecIntegral(t[j + 1], gm, i, false, true)
                            - getRecIntegral(t[j], gm, i, false, true));
                }
            }
            if (isLeg) {
                Betta[i] = 2 * gm * (2 * i + 1) * sum;
            } else {
                Betta[i] = gm * sum;
            }

            System.out.println("Betta [" + i + "]  is - > " + Betta[i]);
            sum = 0;
        }
        return Betta;
    }

    public static double[] findBetta(double t[], double a[],
            double b[], double k, double gm) {
        double[] Betta = new double[(int) (k)];
        double sum = 0;
        for (int i = 0; i < k; i++) {
            for (int j = 0; j < t.length - 1; j++) {
                sum += a[j] * (legRecIntegral(t[j + 1], gm, i)
                        - legRecIntegral(t[j], gm, i))
                        + b[j] * (legRecIntegralT(t[j + 1], gm, i)
                        - legRecIntegralT(t[j], gm, i));
            }
            Betta[i] = 2 * gm * (2 * i + 1) * sum;
            System.out.println("Betta [" + i + "]  is - > " + Betta[i]);
            sum = 0;
        }
        return Betta;
    }

    public static double[] approxFuncOrd(double Betta[], double ti[],
            double gm, double px[]) {
        double[] Kx = new double[ti.length];
        double sum = 0;
        for (int i = 0; i < ti.length; i++) {
            for (double j = 0; j < Betta.length; j++) {
                sum += Betta[(int) j] * leg(ti[i], gm, j);
            }
            Kx[i] = sum;
            System.out.println("Kx[" + i + "] = " + Kx[i] + "   "
                    + " px[" + i + "] = " + px[i]);
            sum = 0;
        }
        return Kx;
    }

    public static double[] approxFuncOrd(double Betta[], double ti[],
            double gm, double px[], boolean isLeg) {
        double[] Kx = new double[ti.length];
        double sum = 0;
        for (int i = 0; i < ti.length; i++) {
            for (double j = 0; j < Betta.length; j++) {
                if (isLeg) {
                    sum += Betta[(int) j] * leg(ti[i], gm, j);

                } else {
                    sum += Betta[(int) j] * lager(ti[i], gm, j);
                }
            }
            Kx[i] = sum;
            System.out.println("Kx[" + i + "] = " + Kx[i] + "   "
                    + " px[" + i + "] = " + px[i]);
            sum = 0;
        }
        return Kx;
    }

    public static double[] findLegPowerSpectrum(double Betta[], double gm,
            double omega, double step) {
        int length = (int) (2 * (omega / step));
        double iter = omega;
        double S[] = new double[length];
        double sum = 0;
        double FiSum = 0;
        omega = 0;
        for (int i = 0; i < S.length; i++) {
            for (int j = 0; j < Betta.length; j++) {
                if (j == 0) {
                    FiSum = 0;
                } else {
                    for (int l = 0; l < j; l++) {
                        FiSum += Math.atan(omega / (2 * gm * (l + 1)));
                    }
                }
                sum += Betta[j] * (1.0 / (2 * j + 1))
                        * Math.cos(Math.atan(omega / (2 * gm * (j + 1))))
                        * Math.cos(Math.atan(omega
                        / (2 * gm * (j + 1))) + 2 * FiSum);
                FiSum = 0;
            }
            S[i] = (1.0 / (Math.PI * gm)) * sum;
            sum = 0;
            omega += step;
        }
        return S;
    }

    public static double[] findLagPowerSpectrum(double Betta[], double gm,
            double omega, double step) {
        int length = (int) (2 * (omega / step));
        double iter = omega;
        double S[] = new double[length];
        double sum = 0;
        double fi = 0;
        omega = 0;
        for (int i = 0; i < S.length; i++) {
            fi = Math.atan((2 * omega) / gm);
            for (int j = 0; j < Betta.length; j++) {
                sum += Betta[i] * (Math.pow(-1, j)) * Math.cos((2 * j + 1) * fi);
            }
            S[i] = (2.0 * Math.cos(fi) / (Math.PI * gm)) * sum;
            System.out.println("S [" + i + "]  - > " + S[i]);
            sum = 0;
            omega += step;
        }
        return S;
    }

    public static double factorial(double val) {
        if (val == 0) {
            return 1;
        } else {
            return val * factorial(val - 1);
        }

    }
}
